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P REFACE  AND  SUMMARY 

A  previous  Memorandum,  RM-5407-PR,  introduced  and  analyzed  a 
simple  multitype  population  model.  That  model  was  special  insofar 
as  the  input  process  was  assumed  Poisson.  This  Memorandum  analyzes 
a  modified  model  that  has  for  its  input  process  an  arbitrary  renewal 
process.  This  modification  allows  far  greater  flexibility  in  fitting 
the  model  to  some  real  situation  than  does  the  case  of  Poisson  inputs. 

S.  C.  Port  is  a  consultant  to  The  RAND  Corporation. 
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A  MULTITYPE  STOCHASTIC  POPULATION  MODEL:  AN  EXTENDED  VERSION 

1 ,  INTRODUCTION 

k 

A  previous  Memorandum  introduced  a  multitype  stochastic  population 

model,  motivated  it  with  examples,  and  analyzed  some  of  its  mathematical 

consequences.  The  purpose  here  is  to  do  the  same  with  a  modified  model 

that  includes  the  original  one  as  a  special  case.  This  new  model  is  the 

same  as  the  previous  one  except  that  the  input  process  is  an  arbitrary 

renewal  process  instead  of  a  Poisson  process.  This  allows  far  greater 

flexibility  in  fitting  the  model  to  a  real  situation. 

The  new  model  is  as  follows.  A  particle  may  be  in  one  of  r  different 

states.  Particles  enter  the  system  in  bunches  that  arrive  according  to 

a  renewal  process  T|( t )  having  waiting  time  f }  ,  which  possesses  a  density 

and  finite  variance.  The  total  number  of  particles  in  the  i-th  bunch  is 

X^,.  The  processes  Tj( t )  and  are  independent.  The  composition  of  the 

i-th  bunch  is  fx.,,  . . .  ,  x .  1;  where  conditional  on  X,,  the  r-vector 
il  ’  ir  i’ 

fx^ . xir^  is  (X^;  ,  ....  or) .  After  arrival  into  the 

system  the  particles  all  act  independently  of  each  other  and  transform 
their  states  according  to  a  semi -Markov  process  having  transition  law 
P, ,(t)  and  final  probabilities,  lim  P.,(t)  -  n,.  We  will  also  assume  that 

iJ  ir  j 


I  |Ptj(t)  *  "jl  dt 


<  a 


As  before,  let  Mj(t)  and  Yj(t)  denote,  respectively,  the  number  of 
particles  that  enter  the  system  in  state  j  by  time  t  and  the  number  of 
particles  in  state  j  at  time  t. 

For  and  negative,  set 


4e  e  3  J  -  VV  V  > 


k 

Sidney  Port,  A  Multitype  Stochastic  Population  Model,  The  RAND 
Corporation,  RM-5407-PR,  September  1967. 
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and  for  o  <  s  <  1 ,  set 

e(s  4)  =  G(s)  , 

and  set  p^t)  -  cYjP^t),  a  ^  1  -  Pj j ( t )  1  -  p2(t)  and  p3(t)  -  o^P^t) 
Then  all  of  the  information  about  the  Y j ( t )  ,  Mj(t)  follows  from  the 
following  integral  equation  that  the  ^2)  satisfy: 

(1.1)  X2)  -  P(WL  >  t)  +  f  g[i  +  (e  1  2  -  l)Pl(t  -  s) 

+  (e  1  -  w*  - s) 

+  -  l)h(l  -  8)]  Vs<V  X2>dP(Wl  *  S>* 

This  equation  may  be  derived  by  the  usual  renewal  type  argument 
and  is  left  to  the  reader.  We  note  that  if  the  are  exponentially 
distributed  with  mean  l/n,  then  the  equation  becomes 

W  V  ’  e"Ut[l  +  *  J  £(•>«,*“"  dal  , 

and  direct  substitution  shows  that 

(1.2)  «x <XX*  X2)  -  exp(p,  f  [G(s)  -  n  ds}  . 

o 

Here 

G(s)  -  g!  1  +  (e  1  2  -  l)Pl(a)  +  (e  1  *  i)p2(8>  +  (e  2  '  !)p3(8)1  • 

Setting  XL  -  19^  X2  -  102  in  (1.2)  and  comparing  with  Eq.  (2.3)  of 
RM-5407-PR,  we  see  that  these  are  the  same.  Thus  in  the  case  of  Poisson 
arrival,  the  Integral  Eq.  (1.1)  becomes  the  basic  equation  of  RM-5407-PR 
We  will  state  and  prove  our  general  results  for  processes  in 
continuous  time.  Naturally,  analogous  results  hold  in  discrete  time. 
Unfortunately,  one  cannot  solve  Eq.  (1.1)  in  the  general  case.  However, 
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we  can  use  it  to  derive  the  behavior  of  the  moments  of  Y^(t)  and  Mj(t) 
Only  the  first  two  moments  will  be  of  interest  to  us.  Set  ^  *  EX^ , 

“  EX^(X^  -  1)  ,  V  “  (EW^)  *■.  Under  our  assumptions,  ET)(t)  -  ty  —  I  ■ 
[Var  W1  -  (EW1)21y2/2>  t  -  ®. 

Theorem  1.  Under  the  assumptions  of  the  model  [Mj(t)/t,  Y^(t)/t] 
converges  <  ^  probability  to  ( ',u  .  Mot  cover  , 

(1.3)  EY j ( t )  -  Mi^jt  +  n^YJ  +  I!  I)  +  o(l)  , 


where 


r 

QQ 

u«lp4J(t>  ■  "Jdt’ 


(1.4)  Var  Y  j(  t)  =  f2u,2Yn2I  +  p^n2y  +  p^H^t  +  o(t)  , 

(1-5)  EM^(t)  «  YUjCy  +  HjOjl  +  °(1), 

and 


(1.6)  Var  M^(t)  -  [2u.2a2Yl  +  u(2)aj^  +  p^Ylt  +  o(l). 


Theorem  2.  Under  the  assumption  of  the  model, 


cov  [Mj( t) ,  y j ( t ) ]  -  [^Y^ry  +  (n(2)  +  upojnjYlt  +  o(t). 


Let  us  now  compare  our  general  results  with  those  for  the  special 
case  when  is  exponential;  i.e.,  T|(t)  is  Poisson  with  the  same 
parameter  values.  Then  the  first  term  in  expansion  of  the  means  is 
the  same  and  the  second  term  differs  by  the  factor  involving  I.  Like¬ 
wise,  in  the  asymptotic  expansion  of  the  variance,  we  see  that  there 
is  an  additional  term  involving  I. 
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The  exponential  waitin  time  is  a  special  case  of  the  gamma.  If 
has  the  gamma  density, 


u  , 

,  .  .  a  u-I  -<yx 
f  (x)  »  ;■  —  x  e  , 

r(u)  ’ 


a  >  o  ,  u  >  o  , 


then 


EW 


1 


u(u  +  1) 
2 
a 


and  thus 


2  2 

„  _  u(u  +1)  „u  u  -  u  _  u(l  -  u) 

1  2  2  “  2  “  2 
a  a  a  o 

Thus  I  >o  if  u  <  1,  I  ”o,  u  =  1,  and  I<o  if  u  >  1.  The  case  u  *  1 
is  the  exponential  distribution. 

An  interesting  special  case  occurs  when  the  arrival  times  are  con¬ 
stant.  Here  it  is  better  to  formulate  the  model  in  discrete  time,  and 
it  is  convenient  to  take  one  time  unit  as  the  time  between  successive 
group  arrivals.  Then  if  Yj(n)  is  the  number  of  particles  in  state  j 
at  time  n,  we  see  that 


V">  ■  x„i +  •• 


+  X 


nn 


where  Xnl  denotes  the  number  of  particles  in  state  j  at  time  n  which 
arrive  at  time  i.  Thus 


Now 


EX 


"J 


Pn‘l 
krkj 


M-jTTj  +  o(l)  , 


< 
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and 


Va 


r  :<n,  .  (Va  r  *,>(]£  Vki')  *  Vkj'X1  '  S  Vkj‘) 


+“lFJ  +°(1)' 


We  therefore  see  that 


Var  Yj(n)  -  nf(Var  ( 1  -  IL)]  +  o(n). 

Thus ,  computing  the  characteristic  function,  we  have 

E  axpp(Yj(n)  -  EY.OD^V.r  Vj(«>j  -  ^  fl  ■  4  ♦  .(£> 

Hence  we  have  established  the  following. 

Theo rem  3 .  For  the  case  of  constant  group  arrivals  the  number 
of  part  la  Is  Y^fn)  in  state  j  at  time  n  normalized  by  its  mean  and 
variance  is  asymptotically  normally  distributed. 
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2,  PROOF  OF  THEOREM  1 

We  first  consider  Y  (t).  Setting  X  =  o  in  (1.1)  and  setting 

El  e  J  1  =  $  (>  )  ,  vie  see  that 
l  t 

if(\)  -  1  -  F(  t)  +  \  g!"  1  +  (e'1  -  1)  p(t  -  s)\t_s(X)  dF(s), 


where 

r 

p(t)  =  pL(t)  +  p3(t>  =  Xj  Q,ipij(t>* 


and  F( t )  =  P(W^  £  t) .  A  differentiation  on  X  yields 

(2.1)  =  f  G’(...)eXP(t  -  b)*  dF(s)  +  f  G( . .  .)*J  dF(s) . 

^  o  *o 

Let  *  F.Xj  and  *  EX^(X1  -  1).  Setting  X  =  o  in  (2.1)  yields 

(2.2)  EY.(t)  =  n  [  p(t  -  s)  dF(s)  +  f  EY . ( t  -  s)  dF(s). 

J  1  Jo  Jo  3 

This  integral  equation  is  readily  solvable  by  Laplace  transforms.  We 
will  illustrate  the  technique  here  and  henceforth  leave  such  details 
to  the  reade".  Denote  Laplace  transforms  by  *  and  let  the  argument  of  the 
transform  function  be  X.  Then  (2.2)  yields 


V**.  A  A  A 

EY j  -  nrpF  +  EY^F 


Hence 


Thus 


(2.3; 


A  A 


EY  (t )  »n,  P(t  -  s)P(W  +  ...  +  W  c  d«) 

j  1  Pi  'o  1 


p(t  -  s)  dN(s)  , 
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where  here  and  in  the  following  N(t)  *  ETl(t)  is  the  renewal  function. 

Under  our  assumptions  on  the  waiting  times  it  is  known  that 
t 

N(t)  =  [*  n(s)  dt  ,  where  n(s)  is  the  renewal  density,  and  that 

rm 

I  |  n(s)  -  1/EW  I  ds  <  *  . 
o 

A  Iso 

«  EW?  -  2(EW  )2  Var  W  -  (EW  )2 

I  •  |  [n(s)  -  1/EW  ]  ds  =  — - - =  - - - x — —  . 

o  (EWp  2(EUl) 

Set  Y  ■  (EWj)  From  (2.3)  we  see  that 

.I- 

EY^(t)  -  j  -  vi-1  |  [p(t  -  s)n(s)  -  Yflj]  ds 

*  U1  j  '  S)  *  n^n(s>  ds 

p* 

+  p.,n  [n(s)  -  yl  ds. 
d  '  o 

Now 

ft 

U.n  \  [n(s)  -  Yl  *  i^H.1  +  °(1),  t  - 

J  o 

Also 


—  w  -V 

1  [p(t  -  s)  -  n  ]n(u)  ds  -  j  [p(t  -  s)  -  IT  IF n(s)  -  y]  ds 
o  J  o  J 

ft 

+  Y  [p(s)  -  n  1  ds. 

J  O  1 

The  first  integral  on  the  right  is  clearly  o(l)  as  t  •*  ®,  and  the 
second  is  YJ  +  o(l),  where 


J 


ITjt  ds. 


\ 


Hence 


rt 

U,  [p(t  -  s)  -  n.]n(s)  ds  =  p,  VJ  +  o(l). 
*•  J  1 


Thus 


(2.4) 


EY 


^(t)  -  =  m^yj  +  n1n^i  +  o 


(i). 


We  must  now  analyze  the  second  moment  EYj(t)  .  Differentiating 
(2.1)  on  \  and  setting  X  =  o  yields  the  equation 

EY  (t)2  =  2  f  n  p(t  -  s ) EY  ( t  -  s)  dF(s)  +u  f  p2(t  -  s)  dF(s) 
^  o  ^  o 

■'t  r1"  2 

+  p1  |  p(t  -  s)  dF(s)  +  I  EY .  (t  -  s)  dF(s)  . 


Solving,  we  find  that 

.2 


(2.5) 


EY 


ds 


.(t)  =  2p  [  p(t  -  s)EY  (t  -  s)n(s) 

J  1  J0  J 

+  W2  I  P  (t  -  s)n(s)  ds  +  EYj(t) 


To  establish  our  results  we  make  good  use  of  the  following. 

Lemma  2.1.  Let  a(t)  and  b(t)  be  nonnegative  functions  such  that 

00 

a(t)  -  aQ ,  b(t)  -  bQ  are  in  L^(o,  ■) ,  and  set  =  J  [a(t)  -  a]  dt , 

=  f  [b(t)  -  b^  dt.  Suppose  M(t)  ^  o,  M(t)  t  ",  t  - 
M(t)  -  tM  -  Ju ,  t  -  ».  Then 

O  M  - 


and 


(2.6) 


pt  a  b  M  t" 

j  b(t  -  s )M( t  -  s)a(s)  ds  -  —  00 
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Proof .  We  can  write 


.t  a  b  M  t 

b(t  -  s)M(  t  -  s)a(s)  ds  -  -  - 


rt 

b(t  -  s)[M(t  -  s)  -  (t  -  s)M  la(S) 

J  _  o 


+  Mq  [  (t  -  s)a(s)b(t  -  s)  ds  - 


M  a  b  t 

OOP 

2 


b(t)[M(t)  -  tM J  ~  bj  t  -  «,  a ( t)  -  a,  t  -  »  . 

O  M 

Thus 

(2.8)  f  b(t  -  s)[M(t  -  s)  -  (t  -  s)M  la(s)  ds  ~  a  b  J t ,  t  -  •. 

**  q  o  o  o  M 

On  the  other  hand, 

(2.9)  M  (t  -  s)a(s)b(t  -  s)  ds  =  M  |  sb(s)a(t  -  s)  ds 

Jo  °  J0 

■  M  [  sb(o)ra(t  -  s)  •  a  ] 

O  J  o 


+  Moao  f  8(b(s)  -  bQ  +  ^f 
o 


M  a  b  t 


Now  integrating  by  parts  we  find  that 


.t  t 

(2.10)  s [b(s )  -  b  ]  ds  »  t  j  F b(s)  -  b  1 

0  O' 


.S 


-  f  ds  f  [b(x)  -  b  ] 


-o  ° 


“  tJb  ‘  tJb  +  o(t)  =  °<0. 


} 
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Also 

f  sb(s)[a(t  -  s)  -  a  ]  ds  “  [  s[b(s)  -  b  lfa(t  -  s)  -  a  1 
..  „  o  J  o  o 

o  o 

+  b  f  s[a(t  -  s)  -  a  1  . 
o  J  o 

o 

The  first  term  on  the  right  is  clearly  o(t),  and  the  second  term  is 
-t  .t 


'o  J  sfa(t  -  s)  -  a^l  ds  =  b^  f  (t  -  s)[a(s)  -  a  ]  ds 
o  0  0  jo  ° 

ft 

Q  J  [a(s)  -  aQl  -  b^  |  s[a(s)  -  aQ]  ds  =  tbQJa  +  o(t). 


Thus 


(2.11)  j*  sb(s)fa(t  -  s)  -  a^]  ds  =  tb^J^  +  o(t). 

“  o 

Substituting  (2.10)  and  (2.11)  into  (2.9)  yields 

t  M  a  b  t2 

(2.12)  M  |  (t  -  s)a(s)b(t  -  s)  ds  =  °  °--- —  +  M  b^t  +  o(t)  . 

o 

Using  this  and  (2.8),  we  find  from  (2.7)  that  (2.6)  holds. 

2 

We  may  now  proceed  to  analyze  EYj(t).  Taking  b(t)  =  p(t)  , 

M(t)  =>  EYj(t),  and  a(t)  =  n(t)  in  the  lemma  and  using  (2.4)  we  find  that 


2nx  J  P(t  -  s)EYj(t  -  s)n(s) 


ds 


2n 


2-2 

rY  n  U 


rr  ri  2  2  1 

— 2  —  1  +  vn^HjYJ  +  M-xn  jl )  t  +  M’^YTTjlt  J  +  o(t)  . 


From  (2.5)  we  then  find  that 

EY2(t)  -  u2v2n2t2  +  [2^1ynj(>i1Yj  +  n  +  2^vn2i]t 


+  (|A2n2y  +  HjYFjJt  +  o ( t ) . 
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On  the  other  hand  f 

[EY  ^(t)  ]2  =  fEVj(t)  -  n)Ynjt][EYu(t)  -  tijYn.tl 

+  2tYHjki1[EYj(t)  -  tl  +  t2ti.^jY2 

=  ^2n2‘/2t2  +  2vnj^1[^1YJ  +  +  o(t). 

Thus 

(2.13)  Var  Y  (t)  =  EY2(t)  -  [EY.(t)]2 

=  [2ii2vn2I  +  +  o ( t )  . 


Turning  our  attention  now  to  M  (t)  we  may  proceed  in  a  similar 

r  XMi(t)i 

manner.  Setting  X  2  “  °  ln  (1.1)  and  letting  $t(\)  »  Ej^e  J  | 
we  find  that 

♦t(X)  -  1  -  F(t)  +  j  cri  +  <e  -  l)aj]«t.6(>)  dF(s). 
Proceeding  as  before  we  find  that 


EMj(t)  -  N(t) 

and 

EMj(t)2  *  2^0^  j  EMj(t  -  s)n(s)  ds  +  ^2o-2N(t)  +  EMj(t) 

-  2u2Oj  j  N(t  -  s)n(s)  ds  +  n2c*2N(t)  +  |±1cTjN(t). 

In  Leona  2.1  set  b(t)  -  1,  M(t)  -  N(t),  and  a(s)  -  n(a). 

Then 

rl  vzt2 

J  N(t  -  s)n(s)  ds  -  —  +  2Ylt  +  o(t). 


-12- 


Thus 

EMj(t)  =  +  2vltl  +  M'2’>  j  Yt  +  na^Vt  +  o(t), 

since 

TEM^t)]2  =  n2Q2N(t)2  =  tx2a^rN(t)  -  Yt]2  +  2vtrN(t)  -  yt]  +  yV) 

2  2  2  2  2  2 
2(1  j Yl  +  (l^Q jY  t  +  o(l). 

Thus 


(2.14)  Var  Mj(t)  =  [  2y.2»2  Yl  +  +  H<>jY"lt  +  o(t). 


We  may  now  easily  establish  the  law  of  large  numbers  for  Y^(t) 
and  Mj(t).  Indeed,  by  (2.13)  we  know  that  Var  Yj(t)  *  0(t) ,  and  thus 


"I  Var  Y  (t) 

Vll  >  6  S  —22 

J  J  e  t 


-  0(p  ■ 


Likewise  (2.14)  shows  that 


IVar  u, 


This  completes  the  proof  of  Theorem  1. 


3.  PROOF  OF  THEOREM  2 


Comput ing 

we  find  that 

(3.1)  Cov  (Mj(t),  Y^(t))  =  |  ^(2)  =  P3(s>)(Pi(s)  +  P2(s>) 

+  UP1(s)  +  nCp^s)  +  P2(s))  EY^t  -  s) 
+  ^(p^s)  =  p3(s))EMj(t  -  s)]n(s)  ds. 


Consider  the  third  term  on  the  right.  Since  p^  +  p2 
becomes 


this  tern 


In  Lemma  2.1  set  b(t)  *  P^t),  EM^t)  -  M(t),and  a(t)  »  n(t).  Then 
in  view  of  Eq.  (1.5)  we  see  that 
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(3.3) 


2  _  2  2 

m.  ,nsv  a.t 


H  P  (s)EM.(t  -  s)n(s)  ds - - * 

1  J  aj  j  2 


=  2ynj^lt 


+  o(t)  . 


Now 


EY^(t)EMj(t)  -  t2^2®^ 

2  2 
=  ^fy^Y1  VJ  +  T^l]t  +  +  o(t). 

Thus  from  (3.1)  -  (3.3)  we  see  that 

Cov  (Yj(t)  ,  Mj(t)  =  [TikjOjVnjl  +  (u(2)  + 

This  establishes  Theorem  2. 


o(t). 
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